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We combine two aspects of magnetic frustration, multiferroicity and emergent quasi-particles 
in spin liquids, by studying magneto-electric monopoles. Spin ice offers to couple these emergent 
topological defects to external fields, and to each other, in unusual ways, making possible to lift the 
degeneracy underpinning the spin liquid and to potentially stabilize novel forms of charge crystals, 
opening the path to a “magnetic crystallography”. In developing the general phase diagram including 
nearest-neighbour coupling, Zeeman energy, electric and magnetic dipolar interactions, we uncover 
the emergence of a bi-layered crystal of singly-charged monopoles, whose stability, remarkably, is 
strengthened by an external [110] magnetic field. Our theory is able to account for the ordering 
process of Tb 2 Ti 2 07 in large field for reasonably small electric energy scales. 


By providing mechanisms for strong magneto-electric 
coupling, frustration has become a key ingredient in mul- 
tiferroics [l|-0. While the search for high-temperature 
multiferroics is appealing for technological application 
such as memory devices [8[, frustration opens a window 
on novel fundamental properties of magnetic matter at 
low temperature where even weak perturbations can play 
an important role. This holds especially for the collective 
behaviour of spin liquids in a wide range of compounds 
from rare-earth oj and copper 10 41211 oxides to organic 
Mott insulators 13| or iridates 1% [l^ • 

In spin ice materials, the constraints imposed by frus¬ 
tration support an extensively degenerate ground state 
where magnetic fluxes are locally conserved [l^. Such 
flux conservation can be described as a divergence-free 
condition, categorizing the spin ice ground state as a 
Coulomb spin liquid by analogy with Maxwell’s electro¬ 
magnetism ITHlQj. where excitations take the form of 
classical magnetic monopoles (Fig. [2jc-tZ) ( 2 ^ . 

In addition to their magnetic properties, it has been 
recently theorized that magnetic monopoles could 
also carry an electric dipole moment [2l| (Fig. [2jc). 
Here we shall investigate the multiple facets of such 
magneto-electric coupling, as an unexplored generic 
ordering process in rare-earth pyrochlores, able to lift 
the degeneracy of spin liquids and to manipulate topo¬ 
logical excitations in frustrated magnets. Our results 
are double. First of all, we give a precise description 
of the mosaic of competing phases in our multiferroic 
spin ice model (Eq.([T])). We show how a ferromagnetic 
double-layer structure of monopoles (DL) is stabilized 
by electric dipolar interactions and even enhanced by a 
magnetic field in the [110] direction. We then use this 
double-layer structure as a signature of multiferroicity in 
rare-earth oxides able to account for recent experiments 
on the spin liquid candidate Tb 2 Ti 2 07 in a field. 


The model - we consider classical Ising spins S 
aligned with their local easy-axes on the pyrochlore lat¬ 
tice supporting electric moments P induced by magneto¬ 


electric coupling [2l| (Fig. [2]), interacting via nearest 
neighbour spin coupling, Zeeman energy, magnetic and 
electric dipolar interactions 
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FIG. 1. When Dm = h = 0, the electric dipoles stabilize a 
monopole double layer (DL - green, see Fig.l^a), in competi¬ 
tion with all in / all out order (AIAO - red), and the Coulomb 
spin liquid (yellow). The circles (crosses) are the transition 
(crossover) temperatures obtained from Monte Carlo simu¬ 
lations. In the hatched regions, even if T = 0 calculations 
confirm the energetic stability of the double-layer structure, 
the first order nature of the transition prevents full thermal- 
ization of the simulations. The solid lines are upper and lower 
mean field estimates of the boundary. See Appendix for de¬ 
tails on simulations and calculations. 
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{b) 2 in - 2 out 
vacuum 


(c) 3 in - 1 out 
single charge 



{d) AIAO 


double charge 


FIG. 2. (a) Electrically induced ground state of our multifer- 
roic spin ice model, made of alternative bi-layers of positive 
(blue) and negative (red) magnetic charges stacked along a 
[001] axis. The a—chains (indicated by thick bonds) carry a 
saturated [110] magnetisation, {b — d) There are three differ¬ 
ent kinds of configurations for a given tetrahedron: 2 in - 2 
out (vacuum of charge forming the Coulomb spin liquid), 3 
in - 1 out (single magnetic charges carrying an electric mo¬ 
ment P whose direction is dictated by the minority (here out¬ 
ward) spin and independent of the sign of the magnetic charge 
because of time-reversal symmetry, forming the double-layer 
structure), and 4 in (double charges forming the AIAO order). 


where i, j and a, /3 are respectively indices for magnetic 
spins on the pyrochlore lattice and electric dipoles on the 
diamond lattice, and Ve = are the respec¬ 

tive nearest neighbour distances. The nearest neighbour 
vector magnetic Si and electric Pa moments have 
unit length. The size of the moments /i and p is included 
in the energy-scale prefactors 

n _ Mo r, _ u _ 3 /o^ 

’”“47rr^’ ^ ^ ^ 


where po and Sq are respectively the vacuum magnetic 
permeability and electric permittivity and H is the ex¬ 
ternal magnetic field. 

The Hamiltonian is studied via classical Monte Carlo 
simulations, using parallel tempering, worm and single¬ 
spin-flip Metropolis algorithm. The dipolar energies have 
been computed with the Ewald summation [22|, [23|, in 
absence of demagnetization factor in order to develop a 
sample-independent theory All spin configurations 

are given in the Appendix. 


Monopole Double Layer - first of all, what happens in 
absence of magnetic interactions, i.e. Dm = h = J = 0 ? 
We find that electric dipoles induce a bi-layer structure 
of single charges with zero polarization and saturated 
magnetization along the [110] axis (Fig. [2ja). Because 


the electric field is even under time reversal, the appari¬ 
tion of such magnetization has to be spontaneous. This 
configuration is unfrustrated at the nearest-neighbour 
level, whose contribution constitutes 96% of the total 
energy. To our knowledge, magneto-electric coupling is 
the first intrinsie mechanism favouring single magnetic 
charges down to zero temperature. 

Loeal ehemieal potential J - in terms of monopoles, 
J plays the role of a chemical potential favoring the 
Coulomb spin liquid for J < 0 (vacuum of charge) and 
the AIAO double-charge crystal for J > 0 [ 23 , [26| , mak¬ 
ing single charges gapped topological excitations in both 
cases. However, the previously observed J = 0 double¬ 
layer structure turns out to be robust over a large range 
of values for J/D^ G [—2.08 : 0.69] (Fig. [T]), raising the 
question on the nature of the mechanism able to stabilize 
such “excitations”. 

For instance when J <0, creating a pair of sin¬ 
gle charges out of the Coulomb spin liquid costs |4J/3| 
while the energy gain is at most —2i^e/3, making such 
monopole-pair creation unfavorable for J/D^ < —1/2. 
An energetically stable cluster of bi-layered monopoles 
thus needs to get bigger and bigger as J decreases in 
order to minimize its surface-over-volume ratio. The 
need for this kind of nucleation process to seed and grow 
a cluster makes the transition first order and prevents 
full thermalization of the simulations in the vicinity of 
the extensively degenerate Coulomb spin liquid (see yel- 
low/green hatched region in Fig. [T]). As a consequence, 
for J —^De^ an experimental cooling down protocol 
would probably fall out-of-equilibrium; once the mag¬ 
net enters the Coulomb spin liquid with a low density 
of monopoles, it will be difficult to nucleate a big enough 
cluster of magnetic charges to crystallize the double-layer 
structure. Such phenomena also exist for J > 0, but to 
a lesser extent because of the low entropy of the AIAO 
ordered phase (see red/green hatched region in Fig. [1]). 

Hence, the magneto-electric opportunity to stabilize 
monopole excitations comes at the cost of large (free) 
energy barriers and multiple metastable states, which 
can naturally account for strong out-of-equilibrium 
effects in pyrochlores. In order to build a comprehensive 
and experimentally relevant picture of the problem, let 
us now include magnetic dipolar interactions, before 
adding a magnetic field able to tune these energy 
barriers, and finally applying our theory to experiments. 

Long range magnetie dipolar interaetions Dm - since 
the electric polarization is coming from magneto-electric 
coupling, it would be improper to neglect the Dm energy 
scale, especially if we keep in mind rare-earth materi¬ 
als with potentially large magnetic moments. In spin 
ice, magnetic dipolar interactions are responsible for the 
effective Coulomb interactions between monopoles [2 q| . 
The property of “projective equivalence” 0 ensures the 
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quasi-degeneracy of the Coulomb spin liquid in presence 
of magnetic dipolar interactions [ 2 ^ , which is only weakly 
lifted in favour of the 2 in - 2 out long range ordered dipo¬ 
lar spin ice (ODSI) state for T Dm [ 27 I. lisj. 

This is why, by favoring 2 in - 2 out configurations, 
opposing the proximity of same-charge monopoles and 
hindering long-range ferromagnetic order, the Dm inter¬ 
action seems very unfavorable to the DL phase. Indeed 
for negative J, the double-layer phase makes way for the 
ODSI at finite Dm (Fig. [3]). Remarkably, the projec¬ 
tive equivalence, valid for magnetic dipolar interactions 
in the Coulomb spin liquid but not electric ones in the 
DL phase, makes the transition temperature one order of 
magnitude smaller from the DL to the ODSI phase. 

However the nearest neighbour contribution of the Dm 
term decreases the monopole chemical potential which 
is four times bigger for double charges than for single 
ones [ 2 ^ . The counter-intuitive consequence is that mag¬ 
netic dipolar interactions favour DL order for positive 
values of J where it was absent at Dm = 0. Because mag¬ 
netic Coulomb interactions in spin ice are four orders of 
magnitude smaller than between bare electric charges at 
the same distance, a relatively low-energy coupling such 
as De is sufficient to counter-balance the repulsion be¬ 
tween neighbouring magnetic charges, opening the path 
for novel crystal structures made of magnetic monopoles. 




FIG. 4. Zero temperature phase diagram in a [110] field h 
obtained from Ewald summation. The arrows show the evo¬ 
lution of the boundaries as h increases. The DL phase is 
further stabilized by h over the AIAO order, but restricted to 
small values of Dm by the apparition of a monolayer phase of 
monopoles (ML, in cyan). 


Because of the intrinsic quasi-degeneracy of the 2 in - 
2 out configurations, the boundary with the DL phase 
only barely shifts as the ODSI order quickly gives way 
to ODSI[iio] upon increasing h. The AIAO phase, on 
the other hand, is suppressed by the [110] field and 
gives way to the DL phase. The electrically induced 
double-layer structure is thus strengthened by the [110] 
magnetic field to larger J and up to Dm ~ 1.66 Dg. For 
Dm > 1.66 De, the effective Coulomb repulsion between 
same-sign monopoles breaks the DL phase in favour of 
a zincblende or monolayer (ML) structure of monopoles 
with saturated magnetisation along the [111] direction. 


FIG. 3. Phase diagram of J/De and DmjD^ obtained from 
Monte Garlo simulations. The vertical axis is the normal¬ 
ized transition temperature into the AIAO (red), double-layer 
(green) or ODSI (orange) phase. See Appendix for details. 

A [110] magnetic field h - in the absence of 
magneto-electric couplings, such a field polarizes half 
of the spins along the a—chains (Fig. [2ja) while mag¬ 
netic dipolar interactions align the remaining spins 
along antiferromagnetically-ordered /3—chains (see Ap¬ 
pendix) . This 2 in - 2 out state is noted ODSI[iio] • 

By ordering the a—chains the field hinders thermal 
fluctuations which improved the thermalization of simu¬ 
lations, even if the presence of modulated phases could 
not be completely ruled out for all boundaries of the four¬ 
dimensional parameter space. The T = 0 phase diagram 
of the simulated phases was then computed by Ewald 
summation (see Appendix and Fig. m). 


Multiferroicity in rare-earth oxides - our results pro¬ 
vide a clear signature of what to look for in experiments 
and remarkably, this double-layer structure has indeed 
been previously observed in Tb 2 Ti 207 under an external 


[110] magnetic field [32 


But is Tb 2 Ti 207 a good candidate for our theory ? We 
believe so for three reasons. Firstly, the Ising anisotropy 
of Tb^+ ions miii and the pinch points observed in po¬ 


larized neutron scattering |36l-l38| are strong indications 
for underlying spin-ice physics. Also, Tb 2 Ti 2 07 pos¬ 
sesses a giant magnetostriction 00, especially along 
the [110] direction [^, and the stability of the low tem¬ 
perature spin liquid phase has been recently ascribed to 
spin-phonon hybridization of the excitations , via dy¬ 
namical Jahn-Teller coupling [i^. Last but not least, 
while sample dependence seems to be an issue in this 
compound, the double-layer structure in a large [110] field 
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scribed as spin liquid [36|, l37|, l45 
antiferromagnetic correlations 


has been confirmed by independent experiments [32|, 
and is thus a robust feature of Tb 2 Ti 207 . It should in¬ 
deed be noted that the nature of the zero- and low-field 
phases of this material remain under debate - possibly 
because of light stuffing/dilution ^ alternatively de- 

ssy with 

5l\ . Such behaviour is 
reminiscent of another rare-earth pyrochlore, Yb 2 Ti 207 : 
while being remarkably well parametrized under a high 
magnetic field (HU, the zero field properties of Yb 2 Ti 207 
noticeably vary between samples ia-Q, also possibly 
due to light stuffing [^. Given the present low-field 
uncertainty, our goal here is to propose an alternative 
scenario for the high field region and to put a new bench¬ 
mark on the 15-year-old puzzle that is Tb 2 Ti 207 . 

Because of a complex single-ion crystal field 113,153- 
13^, the effective size of the Tb^+ magnetic moments is 
not fix ed, but /r = 6 /i^ is a good estimate at low temper¬ 
ature [sj and high field With = 3.59A, Eq. [2] 
gives Dm = 0.48 K. From our theory, the double-layer 
structure can then be stabilized for > D^/ 1.66 ~ 
0.30 K, corresponding to an electric moment of ^ 2.10“^^ 
C.m and a displacement of oxygen ions of ^ 0.6 pm, 
which are reasonable estimates for multiferroics @, 0 . 
Using the parametrization of with J = 2.7 K, our 
multiferroic spin ice model can also explain why the DL 
phase only appears at finite field 33|,[^ (Fig. [5]). In light 
of the sample-dependence issue, it is difficult to push the 
comparison further to low field, where the phase might 
be antiferromagnetic but probably not AIAO [H^, and 
should be separated from the DL structure by a collective 
paramagnet up to 2 Tesla. This is where magnetostric¬ 
tion comes into play as a potential complementary facet 
of our model. 

The [110] direction for a magnetic field has been 
shown to maximize magnetostriction in Tb 2 Ti 207 , 
resulting in field-dependent oxygen displacement |4l| . 
Hence, an external magnetic field stabilizes the DL 
structure not only by suppressing the AIAO order 
(Fig. HI), but possibly also by increasing the electric 
energy scale Dq. In that case, an even wider range of 
parameters {e.g. J = 0.96 K [^) and perturbations 
will lead to a double-layer phase at high field. The 
possibility to include disorder, anisotropic (i^ . 0,0 or 
next-nearest-nei ghb our interactions [63 and quantum 
fluctuations |65l467| to our multiferroic spin ice model 
opens a rich diversity of potential phases to account for 
the yet uncertain phase of Tb 2 Ti 2 07 in low field. In 
particular it is tempting to speculate whether the glassy 
behaviour observed in some samples isl-IHoj might be 
a consequence of the dynamically difficult nucleation 
process discussed in this paper, es pec iall y if Tb 2 Ti 207 


lies close to a spin liquid phase [3£ 


Conclusion - in summary, we have shown how 
magneto-electric coupling can lift the degeneracy of a 
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FIG. 5. Phase diagram parametrized for Tb 2 Ti 2 07 for J = 
2.7 K, Dm = 0.48 K and De = 0.32 K. The [110] magnetic 
field aligns the a chains along the [110] direction, destroying 
the AIAO order in favour for the double-layer structure. All 
error bars are smaller than the dots except for the red/green 
hatched region where simulations were difficult to equilibrate. 
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spin liquid by creating interactions between topological 
excitations. In spin ice these excitations condense into 
a bi-layered monopole crystal strengthened by a [ 110 ] 
magnetic field, a non-trivial example of “magnetic crys¬ 
tallography”. Our theory offers a simple and robust ex¬ 
planation for the ordering of Tb 2 Ti 2 07 in a large [ 110 ] 
field for a reasonably small electric energy scale Dq. 


If we look at the diverse physics emerging from 
itinerant electrons coupled to spin ice (anomalous 
and spontaneous Hall effects in Nd 2 Mo 207 M and 
Pr 2 lr 207 0, non-Kondo resistivity minimum InHii 
and a new kind of quantum criticality in iridates), 
we should expect the coupling to an additional, ferroic, 
degree of freedom to bring a new flavor to spin ice and 
spin liquids, both at equilibrium and dynamically (Toj . 
Experiments on Dy 2 Ti 207 and Ho 2 Ti 207 already 
suggest the presence of magneto-electric effects [77H80|, 
which could be enhanced or even qualitatively modified 
by doping and chemical pressure (sil-ls^ or with the 
inclusion of an electric field. More generally, multiferroic- 
ity offers a promising mechanism to control topological 
defects in magnets. We hope our work will motivate 
further theoretical and experimental investigations of 
multiferroic effects in pyrochlores and spin liquids. 
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cussions. This work was supported by funding from the 
Theory of Quantum Matter unit of the Okinawa Institute 
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APPENDIX 
Ground states energies 

The T = 0 phase diagram of Figs. 1 and 4 in the pa¬ 
per was calculated based on the following energies per 
number of spins N. The spin configurations are given in 
figures [HI [9] and [TOl 


AIAO: E = 

+4.09 Dm 

-J 

( 3 ) 

DL: E = 

0.0455 Dm 

-h/-^ - 0.692 De 

( 4 ) 

ODSI: E = 

-1.95 Dm 

+ t//3 

( 5 ) 

ODSI[iio]: E = 

-1.90 Dm 

-HVq + J/3 

(6) 

ML: E = 

-0.370 Dm 

—hV6 

( 7 ) 

One sees immediately that the DL phase wins 

over 


the ML one for > 0.60D^, and that it is stable for 
J/De G [—2.07; 0.692] when Dm = h = 0. While the 
prefactors for the coupling J and the Zeeman term are 
straightforward to calculate, the remaining terms require 
Ewald summation for a precise estimate. However, many 
of them can actually be calculated analytically to a good 
approximation Hi- Since it provides a useful in¬ 
sight into the analogy between dipoles and monopoles, 
we briefly explain the method in the next section. Please 
note these results are obtained in absence of demagneti¬ 
zation factor. We recover the same values as Yoshida et 
al [33 for the magnetic energies of ODSI and ODSI[iio] 
when including the demagnetization factor of a sphere in 
vacuum. 


Analytical calculation of the energies 

Let Q = riQm denote the magnetic charge on a given 
diamond site where Qm = 2/i/re and n G {—2, —1, 0,1, 2}. 
In presence of magnetic dipolar interactions only, the en¬ 
ergy cost Pn to create a monopole of charge n is (see 
Supplementary Informations of [20|) 



The magnetic interaction between charges is difficult to 
calculate for any random configuration, but if the system 
is charge ordered then it is possible to use the Madelung 
constant of the corresponding crystal structure in order 
to calculate its Coulomb energy Q- Both AIAO and 
ML configurations are ordered in the zincblende struc¬ 
ture with Madelung constant azb = 1-638. With N/2 
diamond sites, the Coulomb energy is 


giving a total energy En = U^- Pn^ 


= 1.53n‘^Dm ( 11 ) 

where n = 1 for the ML and n = 2 for AIAO. This 
value should be compared with a vacuum of charges, i.e. 
the Coulomb spin liquid. Since the degeneracy of the 
Coulomb spin liquid is weakly lifted, the choice for a 
reference energy is somewhat arbitrary. With a lowest 
energy state at — 1.95D^ (ODSI) and the highest energy 
one at — 1.85D^ (fully saturated in the [ 001 ] direction), 
we choose Eref = — 1.90D^ as reference energy and ob¬ 
tain 

AEi = Eml — Eref = 1.53D^vsEi = 1.53D^ 

AE 2 = Eaiao — ^ref = 5.99D^vsE2 = 6.12D^ 

which are in remarkably good agreement. As for the dou¬ 
ble layer structure, we could not find the Madelung con¬ 
stant for such charge ordering. So we calculated it using 
Ewald summation for Coulomb interactions and obtained 
Oi]jL — 0.976. This gives -Edl, Madelung — I- 89 -D 7725 to be 
compared with AEdl = Edl — E^ref = 1-95D^, within 
3% of error. 

Finite temperature simulations 

Now that the zero temperature boundary can be de¬ 
termined exactly from equations dH]) to ([7|), let us turn 
our attention to the finite temperature phase diagram, 
and in particular to the double-layer phase. All simu¬ 
lations were done with parallel tempering, usually with 
1 mK difference between parallel temperatures. To fur¬ 
ther help thermalization, we have developed a variant 
of the worm algorithm for dipolar spin ice [ 2 ^, adapted 
for both electric dipolar interactions and the presence of 
singly charged monopoles. Measurements were made for 
a given number of Monte Carlo steps (MCs) noted tmax^ 
We used two different equilibration processes to com¬ 
pute the error bars: 

• (i) the system is slowly cooled down from high tem¬ 

perature to the temperature T of measurement dur¬ 
ing then it is thermalized at temperature 

T during tmax/^^] for out model, this method pro¬ 
vides a lower bound for the transition temperature. 

• (ii) for any given set of parameters, the ensemble 
of ground states is known exactly (cf. the previous 
sections); the system is then quenched into one of 
the ground state configurations and thermalized at 
temperature T during tmax/^^] this method is bi¬ 
ased towards ordering and offers an upper bound 
for the transition temperature. 
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FIG. 6. Convergence of the transition temperature as a function of measurement time tmax for both equilibration processes 
(i - green) and (ii - blue). For all figures De = 1, h = 0 and N = 432. We confirm that even if the process is very slow, 
simulations converge to the same transition temperature indicated by the dashed line. 
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These two values provide the error bars plotted on Figs. 
1 and 5 of the main text for a system of size N = 1024 
spins. When not visible, the error bars are smaller than 
the symbols. For the 3-dimensional plot of Fig. 3 of the 
main text (no error bars), we only used the equilibration 
process (i) for a system of size N = 432, which provides 
a lower estimate of the stability of the DL phase. For 
Figs. 1, 3 and 5 of the main text, we used tmax = 10^ 
MCS, except in the double layer phase of Fig. 1 where 
tmax = 10^ MCS. 

In Fig. m we show how these two processes converge to 
the same value for four different sets of parameters order¬ 
ing in the double-layer phase (the most difficult ordering 
process in our simulations). Because of the double long- 
range interactions and the very long time of thermaliza- 
tion, big system sizes are difficult to simulate: Fig. [7| 
displays how the transition temperature converges to a 
finite value with increasing linear system size L. 


Husimi tree 

However close to the boundaries, especially with the 
Coulomb spin liquid, if the equilibration process (ii) al¬ 
ways orders in its ground state configuration, the process 

(i) might not be able to find the true ground state and 
will be dominated by the neighbouring phase. This is 
what happens in the hatched regions of Fig. 1 in the 
main text. In that case, we cannot rely solely on simula¬ 
tions to determine the finite temperature phase diagram 
and an analytical approach becomes necessary. 

The out-of-equilibrium region between the DL and 
AIAO is rather narrow, which is why we shall focus on 
the broader one between the Coulomb spin liquid (CSL) 
and the double layer, by estimating the free energy 
of the two phases for Dm = h = 0 (cf. yellow/green 
hatched region in Fig. 1 of the main text). 

Let us first consider the double layer phase. According 
to Eq. m its internal energy is Udl = —0.692De. As 
for the entropy Sdl^ since the transition is strongly 
first order, fluctuations can be neglected in a first 
approximation when compared with the Coulomb spin 
liquid of extensive degeneracy. Thus the free energy is 
Fdl = Udl-TSdl = -0.692 if we fix = 1. 

As for the Coulomb spin liquid, since Dm = h = 0, 
it corresponds to the canonical nearest neighbour spin 
ice model with electric interactions only between singly- 
charged monopoles. To obtain an upper and lower esti¬ 
mate of the boundary between the Coulomb spin liquid 
and the DL phase, we consider the two following cases 

• 1) electric interactions between the dilute 
monopoles are modeled by an effective chemi¬ 


cal potential: the biggest interaction energy gained 
by a pair of monopoles being — 2De/3, we estimate 
an effective chemical potential of —De/3 per 
monopole. 

• 2) electric interactions are neglected: since they 
tend to lower the interacting energies of the dilute 
monopoles, this should give an upper estimate of 
the transition temperature. 


The free energies of both cases can be calculated in the 
Husimi tree approximation 


T 

Fcsl 1 = - ^ log 
T 

Fcsl 2 = - — log 


r4e/5^e/3^e2^j + 3e-2/3J/3- 
2 


( 12 ) 


'4 + e^^'^ + 3e 2^J/3- 
2 


(13) 


Both free energies reproduce the asymptotic limits of 
the Coulomb spin liquid entropy, namely the Pauling 
(T ^ 0+) and paramagnetic (T ^ +00) entropies. 
When compared with Fdl , equations dm and m pro¬ 
vide respectively the lower and upper solid line of Fig. 1 
of the main text for J G [—2.08 : —1.62]. The estimated 
extent of the boundary shift due to electric interactions 
is consistent, within error bars, with our simulations. 
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